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SUMMARY 

A brief description of our improvements and refine- 
ments that led from the CEM95 version of the Cascade- 
Exciton Model (CEM) code to CEM97 and to CEM2k 
is given. The increased accuracy and predictive power 
of the code CEM2k are shown by several examples. To 
describe fission and light- fragment (heavier than 4 He) 
production, the CEM2k code has been merged with 
the GEM2 code of Furihata. We present some results 
on proton-induced fragmentation and fission reactions 
predicted by this extended version of CEM2k. We show 
that merging CEM2k with GEM2 allows us to describe 
many fission and fragmentation reactions in addition to 
the spallation reactions which are already relatively well 
described. Nevertheless, the current version of GEM2 
does not provide a completely satisfactory description 
of complex particle spectra, heavy-fragment emission, 
and spallation yields. We have initiated another ap- 
proach to describe fission, complex particles and frag- 
ment emission by developing further our CEM2k code 
addressing specifically these problems. In this effort, we 
have developed our own universal approximation for in- 
verse cross sections. We have also developed new rou- 
tines to calculate Coulomb barriers and widths of emit- 
ted particles and to simulate their kinetic energy using 
arbitrary approximations for the inverse cross sections. 
To describe fission-fragment production, we have in- 
corporated into CEM2k a thermodynamical model of 
fission by Stepanov. This extended version of CEM2k 
allows us to describe much better complex particle 
emission and many fission fragments, but it is still in- 
complete and needs further work. 

I. INTRODUCTION 

The design of a hybrid reactor system driven with 
a high-current accelerator requires information about 
residual nuclides that are produced by high energy pro- 
tons and secondary neutrons interacting in the tar- 
get and in structural materials. It is both physi- 
cally and economically impossible to measure all nec- 
essary data, which is why reliable models and codes 
are needed. A model with a good predictive power 



both for the spectra of emitted particles and residual 
nuclide yields is the Cascade-Exciton Model (CEM) of 
nuclear reactions 1 . An improved version of the CEM 
is contained in the code CEM95 2 , which is available 
free from the NEA/OECD, Paris. Following an in- 
creased interest in intermediate-energy nuclear data in 
relation to such projects as Accelerator Transmutation 
of nuclear Wastes (ATW), Accelerator Production of 
Tritium (APT), the Spallation Neutron Source (SNS), 
the Rare Isotope Accelerator (RIA), and others, we de- 
veloped a new version of the cascade-exciton model, 
CEM97 3 4 . CEM97 has a number of improved physics 
features, uses better elementary-particle cross sections 
for the cascade model, and due to some significant 
algorithmic improvements is several times faster than 
CEM95. A preliminary version has been incorporated 
into the transport code system MCNPX, although un- 
fortunately the MCNPX team did not include the effi- 
ciency improvements 5 . 

The recent GSI measurements performed using in- 
verse kinematics for interactions of 208 Pb 6 > 7 and 238 U 8 
at 1 GeV/nucleon and 197 Au at 800 MeV/nucleon 9 
with liquid 1 H provide a very rich set of cross sections 
for production of practically all possible isotopes from 
these reactions in a "pure" form, i.e., individual cross 
sections from a specific given bombarding isotope (or 
target isotope, when considering reactions in the usual 
kinematics, p + A). Such cross sections are much eas- 
ier to compare to models than the "camouflaged" data 
from 7-spectrometry measurements. These are often 
obtained only for a natural composition of isotopes in a 
target and are mainly for cumulative production, where 
measured cross sections contain contributions nor only 
from the direct production of a given isotope, but also 
from all its decay-chain precursors. Analysis of these 
new data has motivated us to further improve the CEM 
and to develop a preliminary version of a new code, 
CEM2k 10 ' n , still under development. 

In our original motivation, different versions of the 
CEM codes were developed to reliably describe the 
yields of spallation products and spectra of secondary 
nucleons, without a special emphasis on complex- 
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particle and light-fragment emission or on fission frag- 
ments in reactions with heavy targets. In fact, CEM95, 
CEM97, and the initial version of the CEM2k code sim- 
ulate spallation only and do not calculate the process 
of fission, and do not provide fission fragments and 
a further possible evaporation of particles from them. 
Thus, in simulating the compound stage of a reaction, 
when these codes encounter a fission, they simply tabu- 
late this event (that permits calculation of fission cross 
sections and fissility) and finish the calculation of this 
event without a subsequent treatment of fission frag- 
ments. To be able to describe nuclide production in the 
fission region, these codes have to be extended by in- 
corporating a model of high energy fission (e.g., in the 
transport code MCNPX 5 , where CEM97 and CEM2k 
are used, they are supplemented by Atchison's fission 
model 12 ' 13 , often cited in the literature as the RAL 
model of fission). 

Since many nuclear and astrophysical applications 
require reliable data also on complex particles (gas 
production) and light and/or heavy fragment produc- 
tion, we here address for the first time these questions 
with our improved cascade-exciton model and present 
our preliminary results on the study. We try two dif- 
ferent ways of solving these problems and the results 
from both approaches are presented below after a short 
background on the development of the CEM. 

II. FROM CEM95 TO CEM97 TO CEM2K 
First, we recall the fundamental ingredients of the 
CEM and the main differences between the improved 
cascade-exciton model code CEM97 3,4 and its precur- 
sor, CEM95 2 . The CEM assumes that the reactions oc- 
cur in three stages. The first stage is the IntraNuclear 
Cascade (INC), in which primary particles can be re- 
scattered and produce secondary particles several times 
prior to absorption by, or escape from the nucleus. The 
excited residual nucleus remaining after the emission 
of the cascade particles determines the particle-hole 
configuration that is the starting point for the second, 
pre-equilibrium stage of the reaction. The subsequent 
relaxation of the nuclear excitation is treated in terms 
of the modified exciton model of pre-equilibrium decay 
followed by the equilibrium fission/evaporation stage 
of the reaction. Generally, all three components may 
contribute to experimentally measured particle spectra 
and distributions. 

An important ingredient of the CEM is the crite- 
rion for transition from the intranuclear cascade to 
the pre-equilibrium model. In conventional cascade- 
evaporation models (like ISABEL and Bertini's INC 
used in LAHET 14 , fast particles are traced down to 
some minimal energy, the cutoff energy T cut (or one 
compares the duration of the cascade stage of a re- 
action with a cutoff time, in "time-like" INC models, 
such as the Liege INC 15 ). This cutoff is usually about 



7-10 MeV above the Fermi energy, below which par- 
ticles are considered to be absorbed by the nucleus. 
The CEM uses a different criterion to decide when a 
primary particle is considered to have left the cascade. 

An effective local optical absorptive potential 
W op t.mod.( r ) ' s defined from the local interaction cross 
section of the particle, including Pauli-blocking effects. 
This imaginary potential is compared to one defined by 
a phenomenological global optical model W op t.exp.( r )- 
We characterize the degree of similarity or difference 
of these imaginary potentials by the parameter 

'P | (Wopt.mod. ^^opt.exp. } l^^opt.exp. \ • 

When V increases above an empirically chosen value, 
the particle leaves the cascade, and is then considered 
to be an exciton. Both CEM95 and CEM97 use the 
fixed value V — 0.3. With this value, we find the 
cascade stage of the CEM is generally shorter than 
that in other cascade models. 

The transition from the preequilibrium stage of a re- 
action to the third (evaporation) stage occurs when the 
probability of nuclear transitions changing the num- 
ber of excitons n with An = +2 becomes equal to 
the probability of transitions in the opposite direction, 
with An = —2, i.e., when the exciton model predicts 
an equilibration has been established in the nucleus. 

The improved cascade-exciton model in the code 
CEM97 differs from the CEM95 version by incorporat- 
ing new approximations for the elementary cross sec- 
tions used in the cascade, using more precise values for 
nuclear masses, Q-values, binding and pairing energies, 
using corrected systematics for the level-density param- 
eters, adjusting the cross sections of pion absorption on 
quasi-deuteron pairs inside a nucleus, allowing for nu- 
clear transparency of pions, including the Pauli princi- 
ple in the preequilibrium calculation, and implementing 
significant refinements and improvements in the algo- 
rithms of many subroutines, decreasing the computing 
time by up to a factor of 6 for heavy nuclei, which is 
very important when performing practical simulations 
with transport codes like MCNPX. On the whole, this 
set of improvements leads to a better description of 
the particle spectra and yields of residual nuclei and 
a better agreement with available data for a variety 
of nuclear reactions. Details and examples with some 
results from this work may be found in 3 ' 16 . 

We also make a number of refinements in the calcu- 
lation of the fission channel, as described in 4,17 . Be- 
sides these modifications of the CEM95 code intro- 
duced especially for a better description of fission cross 
sections, we have been further improving the CEM, 
striving for a model capable of predicting different char- 
acteristics of nuclear reactions for arbitrary targets with 
a wide range of incident energies. Modifications made 
for a better description of the preequilibrium, evapo- 
rative, and cascade stages will also affect the fission 
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channel. We have incorporated into the CEM the up- 
dated experimental atomic mass table by Audi and 
Wapstra 18 , the nuclear ground-state masses (where 
data does not exist) and shell corrections by Moller et 
al. 19 , and the pairing energy shifts from Moller, Nix, 
and Kratz 20 into the level-density formula. In addition, 
we have derived a corrected systematics for the level- 
density parameters using the Ignatyuk expression 21 , 
with coefficients fitted to the data analyzed by lljinov et 
al. 22 (we discovered that lljinov et al. used 11/ a/A for 
the pairing energies A (see Eq. (3) in 4 ) in deriving their 
level-density systematics instead of the value of 12/y/A 
stated in 22 ). We also derived additional semiempirical 
level-density-parameter systematics using the Moller et 
al. 19 ground-state microscopic corrections, both with 
and without the Moller, Nix, and Kratz 20 pairing gaps, 
and introduced into the CEM a new empirical relation 
to take into account the excitation-energy dependence 
of the ground-state shell correction in the calculation 
of fission barriers (see 4 ). 

As mentioned in the Introduction, analysis of the re- 
cent GSI measurements 6-9 has motivated us to further 
improve the CEM. The authors of the GSI measure- 
ments performed a comparison of their data to several 
codes, including LAHET 14 , YIELDX 23 , ISABLA (IS- 
ABEL INC code from LAHET followed by the ABLA 24 
evaporation code), CASCADO 25 , and the Liege INC by 
Cugnon 15 , and encountered serious problems; none of 
these codes were able to accurately describe all their 
measurements. Most of the calculated distributions 
of isotopes produced as a function of neutron number 
were shifted toward larger masses as compared to the 
data. While in some disagreement with the measure- 
ments, the Liege INC and the CASCADO codes pro- 
vide a better agreement with the data than LAHET, 
ISABLA, and YIELDX do. Being aware of this situa- 
tion with the GSI data, we decided to consider them 
ourselves, leading to the development of CEM2k. 

First, we calculated the 208 Pb GSI reaction 6 with 
the standard versions of CEM95 and CEM97 and 
determined 10 that though CEM95 describes quite well 
production of several heavy isotopes near the target 
(we calculate p + 208 Pb; therefore 208 Pb is a target, 
not a projectile as in the GSI measurements), it does 
not reproduce correctly the cross sections for lighter 
isotopes in the deep spallation region. The disagree- 
ment increases with increasing distance from the tar- 
get, and all calculated curves are shifted to the heavy- 
mass direction, just as was obtained by the authors 
of the GSI measurements with the codes they tried. 
The results of the CEM97 code are very similar to 
those of CEM95 (see a figure with CEM97 results 
in 26 ). Later on, we performed an extensive set of cal- 
culations of the same data with several more codes 
(HETC 27 , LAHET 14 with both ISABEL and Bertini 
options, CASCADE 28 , CASCADE/INPE 29 , INUCL 30 , 



and YIELDX 23 ) and got very similar results: 26 all codes 
disagree with the data in the deep spallation region, 
the disagreement increases as one moves away from 
the target, and all calculated curves are shifted in the 
heavy-mass direction. 

This means that for a given final element (Z), all 
models predict emission of too few neutrons. Most 
of the neutrons are emitted at the final (evaporation) 
stage of a reaction. One way to increase the number 
of emitted neutrons would be to increase the evapo- 
rative part of a reaction. In the CEM, this might be 
done in two different ways: the first would be to have 
a shorter preequilibrium stage, so that more excitation 
energy remains available for the following evaporation; 
the second would be to have a longer cascade stage, 
so that after the cascade, less exciton energy is avail- 
able for the preequilibrium stage, so fewer energetic 
preequilibrium particles are emitted, leaving more exci- 
tation energy for the evaporative stage. 

One easy way to shorten the preequilibrium stage of 
a reaction in CEM is to arbitrarily allow only transi- 
tions that increase the number of excitons, An = +2, 
i.e., only allow the evolution of a nucleus toward the 
compound nucleus. In this case, the time of the equi- 
libration will be shorter and fewer preequilibrium par- 
ticles will be emitted, leaving more excitation energy 
for the evaporation. Such an approach is used by some 
other exciton models, for instance, by the Multistage 
Preequilibrium Model (MPM) used in LAHET 14 . Cal- 
culations using this modification to CEM97 (see Fig. 
2d in 10 ) provide a shift of the calculated curves in the 
right direction, but only very slightly improve agree- 
ment with the GSI data. 

A second method of increasing evaporation is to en- 
large the cascade part of a reaction; we may either 
enlarge the parameter V or remove it completely and 
resort to a cutoff energy T cut , as is done in other INC 
models. Calculations have shown that a reasonable in- 
crease of V doesn't solve the problem. However, if 
we do not use V, but instead use a cutoff energy of 
T cu t = 1 MeV for incident energies above the pion 
production threshold, the code agrees with the GSI 
data significantly better (see Fig. 2e in 10 ). Using both 
these conditions leads to results that describe the p 
+ 208 Pb GSI data very well. We call this approach 
CEM2k. An example of CEM2k results for the yield of 
Tm, Ir, and Tl isotopes from p + 208 Pb interactions 
compared to the GSI 6 and ITEP 31 data and with pre- 
dictions by CEM95, LAHET-ISABEL, LAHET-Bertini, 
CASCADE, CASCADE/INPE, INUCL, and YIELDX is 
shown in Fig. 1. Similar comparisons for more iso- 
topes may be found in 10 < 26 < 31 . We find that CEM2k 
agrees best with these GSI (and ITEP) data of the 

codes tested here and i n 10 ' 26 ' 31 . 

Finding a good agreement of CEM2k with the iso- 
tope production, we wish to see how well it de- 
scribes spectra of secondary particles in comparison 
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Isotopic distributions of the products in Pb+1 GeV protons: GSI+ITEP+Codes 




Figure 1. Mass distributions for independent pro- 
duction of Tm, Ir, and Tl isotopes from 1 GeV 
protons colliding with 208 Pb. Stars are GSI 
measurements 6 , while squares show ITEP data 31 . 
The calculations are identified as: CEM2k — our re- 
sults, CEM95 2 , LAHET-ISABEL 14 , LAHET-Bertini 14 , 
CASCADE/INPE 29 , CASCADE 28 , INUCL 30 , and 
YIELDX 23 . 

with CEM97. Figure 2 shows examples of neutron 
spectra from interactions of protons with the same 
target, 208 Pb at 0.8 and 1.5 GeV (we do not know 
of measurements of spectra at 1 GeV, the energy of 
the isotope-production data). We see that CEM2k de- 
scribes spectra of secondary neutrons comparably to 
CEM97, even possibly a little better at larger angles. 

So this preliminary version of CEM2k describes both 
the GSI data from 208 Pb interactions with p at 1 
GeV/nucleon and the spectra of emitted neutrons from 
p+ 208 Pb at 0.8 and 1.5 GeV better than its precursor 
CEM97. 

We use CEM2k as fixed from our analysis of the 
208 Pb data 6 ' 7 without further modifications to calcu- 
late the 197 Au 9 and 238 U 8 GSI measurements. An ex- 
ample of the yield of several isotopes from 197 Au cal- 
culated by CEM2k is shown in Fig. 3 together with 
standard CEM97 predictions and calculations by the 
LAHET-Bertini and YIELDX codes from 9 . We see that 
just as in the case of the 208 Pb data, CEM2k agrees 
best with the 197 Au data in the spallation region com- 
pared to the other codes tested here. Several more 
results for 197 Au and 238 U and their detailed discus- 
sion may be found in 10 . 

Besides the changes to CEM97 mentioned above, 
we also made a number of other improvements and 
refinements, such as imposing momentum-energy con- 
servation for each simulated event (the Monte Carlo 
algorithm previously used in CEM provides momentum- 
energy conservation only statistically, on the average, 
but not exactly for the cascade stage of each event); 
using real binding energies for nucleons at the cascade 
stage instead of the approximation of a constant sepa- 
ration energy of 7 MeV used in previous versions of 
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Figure 2. Comparison of measured 32 double differen- 
tial cross sections of neutrons from 0.8 and 1.5 GeV 
protons on Pb to CEM2k and CEM97 calculations. 
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Figure 3. Isotopic distributions of spallation products 
from the reaction 197 Au + p at 800 A MeV from 
mercury to hafnium. Open circles are the GSI data 9 , 
CEM2k (thick solid curves) and CEM97 (thick dashed 
curves) are our present calculations, LAHET-Bertini 
(thin solid curves) and YIELDX (thin dashed curves) 
are calculations from 9 . 
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the CEM; and using reduced masses of particles in 
the calculation of their emission widths instead of us- 
ing the approximation of no recoil used previously. In 
Ref 11 . we have shown that these refinements (which 
involve no further parameter fitting), while only im- 
proving slightly the agreement with the GSI measure- 
ments (and other data on medium and heavy targets), 
are crucial for light targets, especially when calculat- 
ing 4 He and other fragment emission from light nuclei. 
This is especially important for applications where gas 
production is calculated. 

Another improvement important for applications is 
better representing the total reaction cross section. 
Previous versions of CEM (just like many other IIMC- 
type models) calculate the total reaction cross section, 
(Ji n using the geometrical cross section, a geom , and 
the number of inelastic, iVj„, and elastic, N e i, simu- 
lated events, namely: a in = cr geom N m /(N m + N el ). 
This approach provides a good agreement with avail- 
able data at incident energies above about 100 MeV, 
but is not reliable at lower bombarding energies. To ad- 
dress this problem, we have incorporated into CEM2k 
the NASA systematics by Tripathi et al. 33 for all in- 
cident protons and neutrons with energies above the 
maximum in the NASA reaction cross sections, and 
the Kalbach systematics 34 for neutrons of lower energy. 
As shown in Fig. 4, we can describe much better with 
CEM2k the total reaction cross sections (and corre- 
spondingly any other partial cross sections) for n- and 
p-induced reactions, especially at energies below about 
100 MeV. 

III. MERGING CEM2K WITH GEM2 
As a first attempt to describe with CEM2k both 
emission of intermediate-mass fragments heavier than 
4 He and production of heavy fragments from fission, 
we merged CEM2k with the Generalized Evaporation 
Model (GEM) code by Furihata 37 ' 38 . GEM is an 
extension by Furihata of the Dostrovsky evaporation 
model 39 as implemented in LAHET 14 to include up to 
66 types of particles and fragments that can be evap- 
orated from an excited compound nucleus plus a mod- 
ification of the version of Atchison's fission model 12,13 
used in LAHET. Many of the parameters were adjusted 
for a better description of fission reactions when using 
it in conjunction with the extended evaporation model. 
We merged GEM2 (the last update of the GEM code) 
with CEM2k as follows: we calculate the cascade and 
preequilibrium stages of a reaction with our CEM2k, 
then we describe the subsequent evaporation of par- 
ticles and fragments and fission from the remaining 
excited compound nuclei using GEM2. To understand 
the role of preequilibrium particle emission, we per- 
formed calculations of all the reactions we tested both 
with emission of preequilibrium particles and without 
them, i.e., going directly to GEM2 after the intranu- 
clear cascade stage of a reaction described by CEM2k. 
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p-induced reactions on Al calculated by CEM2k and 
CEM97 compared with experimental data compiled by 
Barashenkov 35 and calculations from the HMS-ALICE 

code 36 . 
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A very detailed description of the GEM, together 
with a large amount of results obtained for many re- 
actions using GEM coupled either with the Bertini or 
ISABEL INC models in LAHET may be found in 37 ' 38 . 
Therefore, we present here only the main features of 
GEM, following mainly 38 and using as well useful in- 
formation obtained in private communications from Dr. 
Furihata. 

A. Evaporation Model in GEM 
Furihata did not change in GEM the general algo- 
rithms used in LAHET to simulate evaporation and 
fission. The decay widths of evaporated particles and 
fragments are estimated using the classical Weisskopf- 
Ewing statistical model 40 . In this approach, the decay 
probability Pj for the emission of a particle j from a 
parent compound nucleus i with the total kinetic en- 
ergy in the center-of-mass system between e and e + de 
is 

Pj{t)de = g a lnv {e) — ede, (1) 

where E [MeV] is the excitation energy of the parent 
nucleus i with mass and charge Z it and d denotes 
a daughter nucleus with mass Ad and charge Zd pro- 
duced after the emission of ejectile j with mass Aj and 
charge Zj in its ground state. a inv is the cross section 
for the inverse reaction, pi and pd are the level densi- 
ties [MeV] -1 of the parent and the daughter nucleus, 
respectively, gj = (2Sj + \)nij /tt 2 Ti 2 , where Sj is the 
spin and mj is the mass of the emitted particle j. The 
Q-value is calculated using the excess mass M(A,Z) 
as Q = M(Aj, Zj)+M(A d , Z d )-M{A h ZA. In GEM, 
four mass tables are used to calculate Q-values, ac- 
cording to the following priority: (1) the Audi-Wapstra 
mass table 18 , (2) theoretical masses calculated by 
Moller et al. , (3) theoretical masses calculated by 
Comay et al. 41 , (4) the mass excess calculated using 
the old Cameron formula 42 . As does LAHET, GEM 
uses Dostrovsky's formula 39 to calculate the inverse 
cross section a inv for all emitted particles and frag- 
ments 

(?inv{e) = VgOL (l + ^j , (2) 

which is often written as 

( \ — / u a Cn ^y + k/ e ) f° r neu trons 
VinvW - | agC .^ _ yj^ for charged part i c | e s , 

where a g = ttR 2 [fm 2 ] is the geometrical cross section, 
and 

V = k J Z J Z d e 2 /R c (3) 
is the Coulomb barrier in MeV. 



The new ingredient in GEM in comparison with LA- 
HET which considers evaporation of only 6 particles 
(n, p, d, t, 3 He, and 4 He) is that Furihata included the 
possibility of evaporation of up to 66 types of particles 
and fragments and incorporated into GEM several sets 
of parameters b, Cj, kj, Rb, and R c for each particle. 

The 66 ejectiles considered by GEM for evaporation 
are selected to satisfy the following criteria: (1) iso- 
topes with Zj < 12; (2) naturally existing isotopes or 
isotopes near the stability line; (3) isotopes with half- 
lives longer than 1 ms. All the 66 ejectiles considered 
by GEM are shown in Table 1. 

Table 1. The ejectiles considered by GEM 



Z 3 Ejectiles 






n 














1 


P 


d 


t 










2 


3 He 


4 He 


6 He 


8 He 








3 


6 Li 


7 Li 


8 Li 


9 Li 








4 


7 Be 


9 Be 


10 Be 


n Be 


12 Be 






5 


8 B 


10 B 


n B 


12 B 


13 B 






6 


10 C 


n C 


12 C 


13 C 


14 Q 


15 C 


16 C 


7 


12 N 


13 N 


14 N 


15 N 


16 N 


17 N 




8 


14Q 


15Q 


16Q 


17 


18Q 


19Q 


20 


9 


17 F 


18 F 


19 F 


20 f 


21 F 






10 


18 Ne 


19 Ne 


20 Ne 


21 Ne 


22 Ne 


23 Ne 


24 Ne 


11 


21 Na 


22 Na 


23 Na 


24 Na 


25 Na 






12 


22 Mg 


23 Mg 


24 Mg 


25 Mg 


26 Mg 


27 Mg 


28 Mg 



GEM includes several options for the parameter set 
in expressions (2,3): 

1) The "simple" parameter set is given as c n = Cj — 

k 3 = 1 b = 0, and R b = R c = MA 1 / 3 + A d /3 ) [fm]; 
users need to input r . 

2) The "precise" parameter set is used in GEM as 
default, and we use this set in our present work. 

A) For all light ejectiles up to a (Aj < 4), the pa- 
rameters determined by Dostrovsky et al. 39 are used in 
GEM, namely: c„ = 0.76 + c a A~ 1/3 , b = {b a A~ 2/S - 

0.050)/(0.76 + c a A^ 1/3 ) (and b = for A d > 192), 
where c a = 1.93 and b a = 1.66, c p = 1 + c, 
c d = 1 + c/2, c t = 1 + c/3, C3 He = c a = 0, k p = k, 
k d = k + 0.06, k t = k + 0.12, k?, He = k a - 0.06, 
where c, k, and k a are listed in Table 2 for a set of 
Z d - Between the Z d values listed in Table 2, c, k, and 
k a are interpolated linearly.The nuclear distances are 
given by Rb = 1.5A 1 / 3 for neutrons and protons, and 
1.5{A d /s + A 1 / 3 ) for d, t, 3 He, and a. 

Table 2. k, k a , and c parameters used in GEM 



z d 


k 




c 


< 20 


0.51 


0.81 


0.0 


30 


0.60 


0.85 


-0.06 


40 


0.66 


0.89 


-0.10 


> 50 


0.68 


0.93 


-0.10 
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The nuclear distance for the Coulomb barrier is ex- 
pressed as R c = R d + Rj, where R d = r^A 1 / 3 , 
Tq = 1.7, and Rj = for neutrons and protons, and 
Rj = 1.2 for d, t, 3 He, and 4 He. We note that several 
of these parameters are similar to the original values 
published by Dostrovsky et al. 39 (and used, for example 
in our CEM codes) but not exactly the same. Dostro- 
vsky et al. 39 had c a = 2.2, b a = 2.12, and rg = 1.5. 
Also, for the k, k a , and c parameters shown in Table 
2, they had slightly different values, shown in Table 3. 



Table 3. k p , c p , k a , and c a parameters from Ref. 



:!') 



z d 


kp 


Cp 






10 


0.42 


0.50 


0.68 


0.10 


20 


0.58 


0.28 


0.82 


0.10 


30 


0.68 


0.20 


0.91 


0.10 


50 


0.77 


0.15 


0.97 


0.08 


> 70 


0.80 


0.10 


0.98 


0.06 



B) For fragments heavier that a (Aj > 4), "the 
precise" parameters of GEM use values by Matsuse et 
al. 43 , namely: c 3 ■ = k = 1, R b = R Q (Aj) + R (A d ) + 
2.85 [fm], R c = R (Aj) + R (A d ) + 3.75 [fm], where 
R (A) = 1.12AV3 - 0.86A- 1 / 3 . 

3) The upgraded version of GEM realized in the code 
GEM2 contains two other options for the parameters 
of the inverse cross sections. 

A) A set of parameters due to Furihata for light 
ejectiles in combination with Matsuse's parameters for 
fragments heavier than a. Furihata and Nakamura de- 
termined kj for p, d, t, 3 He, and a as follows: 44 

kj = ci \og(Z d ) + c 2 \og(A d ) + c 3 . 

The coefficients c\, c 2 , and c 3 for each ejectile are 
shown in Table 4. 

Table 4. ci, c 2 , and c 3 for p, d, t, 3 He, and a from 44 



Ejectile 


Cl 


C2 


C3 


P 


0.0615 


0.0167 


0.3227 


d 


0.0556 


0.0135 


0.4067 


t 


0.0530 


0.0134 


0.4374 


3 He 


0.0484 


0.0122 


0.4938 


a 


0.0468 


0.0122 


0.5120 



When these parameters are chosen in GEM2, the 
following nuclear radius R is used in the calculation of 



V and (j g : 



R 





1.2 
2.02 
2.42 
2.83 
3.25 
k 1.414AV 3 + 1 



for A = 1 , 
for 2 < A < 4 
for 5 < A < 6 
for A = 7 , 
for A = 8 , 
for A = 9 , 
for A > 10 . 



B) The second new option in GEM2 is to use Fu- 
rihata's parameters for light ejectiles up to a and the 
Botvina et al. 45 parameterization for inverse cross sec- 
tions for heavier ejectiles. Botvina et al. 45 found that 
<r inv can be expressed as 



G 9 | 



(1-V/e) 

V+l 



fore> V + l [MeV], 
fore<y+l [MeV], 



(4) 



where 



a = 0.869 + 9M/Zj, 



V = 



ZjZ d 



2.173 



1 



](Af+Ary 



6.103 x 10- 3 Z ? Z d rr . 

[fm] 



1 + 9.443 x 10- 3 Z 7 Z d 



The expression of <Ji nv for e < V + 1 shows the 
fusion reaction in the sub-barrier region. When using 
Eq. (4) instead of Eq. (2), the total decay width for a 
fragment emission can not be calculated analytically. 
Therefore, the total decay width must be calculated 
numerically and takes much CPU time. 

The total decay width Tj is calculated by integrating 
Eq. (1) with respect to the total kinetic energy e from 
the Coulomb barrier V up to the maximum possible 
value, (E — Q). The good feature of Dostrovsky's 
approximation for the inverse cross sections, Eq. (2), 
is its simple energy dependence that allows the analytic 
integration of Eq. (1). By using Eq. (2) for cn nv , the 
total decay width for the particle emission is 

The level density p(E) is calculated in GEM according 
to the Fermi-gas model using the expression 46 



P(E) 



it exp(2y/a(E - 5)) 
12 a 1 / 4 (E-S) 5 / 4 ' 



(6) 



where a is the level density parameter and S is the 
pairing energy in MeV. As does LAHET, GEM uses 
the 5 values evaluated by Cook et al. 47 For those val- 
ues not evaluated by Cook et al., S's from Gilbert 
and Cameron 46 are used instead. The simplest option 
for the level density parameter in GEM is a = A d /8 
[MeV -1 ], but the default is the Gilbert-Cameron-Cook- 
Ignatyuk (GCCI) parameterization from LAHET: 14 



1 



a = a- 



ai \ 1 



1 



(7) 
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where u = 0.05(E — S), and 

ai = (0.1375 - 8.36 x lCT 5 ^) x A d 



a = 



A d /8 

A d (a' + 0.009175") 



for Z d < 9 or N d < 9, 
for others. 



For deformed nuclei with 54 < Z d < 78, 86 < Z d < 
98, 86 < N d < 122, or 130 < N d < 150, a' = 0.12 
while a' = 0.142 for other nuclei. The shell correc- 
tions S is expressed as a sum of separate contributions 
from neutrons and protons, i.e. S — S(Z d ) + S(N d ) 
from 46 ' 47 and are tabulated in 38 . 

The level density is calculated using Eq. (6) only for 
high excitation energies, E > E x , where E x = U x + 8 
and U x = 2.5 + 150/ A d (all energies are in MeV). At 
lower excitation energies, the following 46 is used for the 
level density: 

P( E ) = ^cM(E-E )/T), (8) 

where T in the nuclear temperature defined as 1/T = 
\faJTT x — 1.5/U X . To provide a smooth connection of 
Eqs. (6) and (8) at E = E x , E is defined as Eq = 
E x - T(\ogT - 0.25 loga - 1.25 log £4 + 2*JaU x ). 

For E - Q - V < E x , substituting Eq. (8) into 
Eq. (5) we can calculate the integral analytically, if we 
neglect the dependence of the level density parameter 
a on E: 

irgjO-ga 



Tj l2 Pi {E) 



{h(t,t) + (J3 + V)I (t)}, 



(9) 



where Io(t) and h(t,t x ) are expressed as 



-E /T 



1) 



I {t) = e 

h(t,t x ) = e- E °/ T T{(t-t x + l)e t * -t-l}, 

where t = (E - Q - V)/T and t x = E x /T. For 
E — Q - V > E x , the integral of Eq. (5) cannot be 
solved analytically because of the denominator in Eq. 
(6). However, it is approximated as 

T] = nME) [h{t > tx) + /3(s ' 3xr + ^ + v ^ 

x {Io(t x )-h(s,s x )e s }], (10) 
where l2(s,s x ) and l3(s,s x ) are given by 

I 2 (s, s x ) = 2V2{ S - 3/2 + 1.5s- 5/2 + 3.75s- 7/2 

-( S - 3 / 2 + 1.5 S - 5 /2 + 3. 75s -7/ 2)eS ,- s}; 



h(a, s x ) - (V^)- 1 ^- 172 + 4s- 3 / 2 + 13.5s- 5 / 2 
+ 60.0s" 7/2 + 325.125s" 9/2 



{( s ~ S x) S x 



-3/2 



0.25<)s- 



(1.5s 2 

7/2 



0.5s 2 x )s x 



5/2 



+ (3.75s 2 
+ (12.875s 2 + 0.625s 2 ) S:i : 9/2 
+ (59.0625s 2 + 0.9375s 2 )s 2 : 11/2 
+ (324.8s 2 + 3.28s 2 )s- 13 / 2 }e^" 



with s = 2y / o(E — Q — V — S) and s x — 

2y/a(E x -6). 

The ejectile j to be evaporated is selected in GEM by 
the Monte Carlo method according to the probability 
distribution calculated as Pj = Tj/Y^j^j, where Tj 
is given by Eqs. (9) or (10). The total kinetic energy 
e of the emitted particle j and the recoil energy of the 
daughter nucleus is chosen according to the probability 
distribution given by Eq. (1). The angular distribution 
of ejectiles is simulated to be isotropic in the center- 
of-mass system. 

According to Friedman and Lynch 48 , it is important 
to include excited states in the particle emitted via the 
evaporation process along with evaporation of parti- 
cles in their ground states, because it greatly enhances 
the yield of heavy particles. Taking this into consider- 
ation, GEM includes evaporation of complex particles 
and light fragments both in the ground states and ex- 
cited states. An excited state of a fragment is included 
in calculations if its half-lifetime 7\/ 2 (s) satisfies the 
following condition: 



Tl/2 

~W2 



> 



(11) 



where V* is the decay width of the excited particle (res- 
onance). GEM calculates T* in the same manner as 
for a ground-state particle emission. The Q-value for 
the resonance emission is expressed as Q* = Q + E* , 
where Ej is the excitation energy of the resonance. 
The spin state of the resonance S* is used in the cal- 
culation of gj, instead of the spin of the ground state 
Sj. GEM uses the ground state masses mj for ex- 
cited states because the difference between the masses 
is negligible. 

Instead of treating a resonance as an independent 
particle, GEM simply enhances the decay width Tj of 
the ground state particle emission as follows: 



(12) 



where T° is the decay width of the ground state particle 
emission, and T™ is that of the nth excited state of the 
particle j emission which satisfies Eq. (11). 

The total kinetic energy distribution of the excited 
particles is assumed to be the same as that of the 
ground state particle emission. S* , E* , and Ti/ 2 
used in GEM are extracted from the Evaluated Nuclear 
Structure Data File (ENSDF) database maintained by 
the National Nuclear Data Center at Brookhaven Na- 
tional Laboratory 49 . 

Note that when including evaporation of up to 66 
particles in GEM, its running time increases many times 
compared to the case when evaporating only 6 parti- 
cles, up to 4 He. The major particles emitted from an 
excited nucleus are n, p, d, t, 3 He, and 4 He. For most 
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of the cases, the total emission probability of particles 
heavier than a is negligible compared to those for the 
emission of light ejectiles. To keep the GEM running 
time to a reasonable level, calculations of probability 
of emission of particles heavier than a are done only 
for 5% of simulated evaporations for Ai > 40, 7% of 
simulated evaporations for 30 < Ai < 40, and 30% of 
simulated evaporations for 20 < Ai < 30. This empir- 
ical criterion was determined by Furihata from actual 
simulations 38 . 

B. Fission Model in GEM 

The fission model used in GEM is based on Atchison's 
model 12,13 as implemented in LAHET 14 , often referred 
in the literature as the Rutherford Appleton Labora- 
tory (RAL) model, which is where Atchison developed 
it. There are two choices of parameters for the fis- 
sion model: one of them is the original parameter set 
by Atchison 12 ' 13 as implemented in LAHET 14 , and the 
other is a parameter set evaluated by Furihata 37 ' 38 . 

B.l. Fission Probability. The Atchison fission 
model is designed to only describe fission of nuclei with 
Z > 70. It assumes that fission competes only with 
neutron emission, i.e., from the widths Tj of n, p, d, 
t, 3 He, and 4 He, the RAL code calculates the proba- 
bility of evaporation of any particle. When a charged 
particle is selected to be evaporated, no fission compe- 
tition is taken into account. When a neutron is selected 
to be evaporated, the code does not actually simulate 
its evaporation, instead it considers that fission may 
compete, and chooses either fission or evaporation of 
a neutron according to the fission probability Pf. This 
quantity is treated by the RAL code differently for the 
elements above and below Z = 89. The reasons Atchi- 
son split the calculation of the fission probability Pf 
are: (1) there is very little experimental information on 
fission in the region Z = 85 to 88, (2) the marked rise 
in the fission barrier for nuclei with Z 2 /A below about 
34 (see Fig. 2 in 13 ) together with the disappearance 
of asymmetric mass splitting, indicates that a change 
in the character of the fission process occurs. If ex- 
perimental information were available, a split between 
regions about Z 2 /A s=a 34 would more sensible 13 . 

1) 70 < Zj < 88. For fissioning nuclei with 70 < 
Zj < 88, GEM uses the original Atchison calculation 
of the neutron emission width T n and fission width Tf 
to estimate the fission probability as 



Pf = 



f 



1 



r / + r„ i + r„/ry 



(13) 



Atchison uses 12 ' 13 the Weisskopf and Ewing sta- 
tistical model 40 with an energy-independent pre- 
exponential factor for the level density (see Eq. (6)) 
and Dostrovsky's 39 inverse cross section for neutrons 
and estimates the neutron width T n as 



r„ = 0.352(1. 68J + IMA 1 / 3 J x 
+A 4 2/3 (0.76Ji-0.05J )), 



(14) 



where Jo and J\ are functions of the level density pa- 
rameter a n and s n (= 2y / a n (E — Q n — 5)) as 



Jo 



(gn - l)e«" + 1 
2a n 



Ji = 



(2s 2 - 6s n + 6)e s " +a 2 -6 
8a 2 



Note that the RAL model uses a fixed value for the 
level density parameter a n , namely 



a n = (A, - l)/8, 



(15) 



and this approximation is kept in GEM when calculat- 
ing the fission probability according to Eq. (13), though 
it differs from the GCCI parameterization (7) used in 
GEM to calculate particle evaporation widths. The fis- 
sion width for nuclei with 70 < Zj < 88 is calculated 
in the RAL model and in GEM as 



r/ = 



(Sf - l)e*/ + 1 
a f 



(16) 



where Sf = 2y/a,f{E — Bf — S) and the level density 
parameter in the fission mode a/ is fitted by Atchison 
to describe the measured Tf/T n as: 13 



«/ 



1.08926 + 0.01098(x - 31.08551) 2 ) , (17) 



and x = Z 2 /A. The fission barriers B f [MeV] are 
estimated as 



Z 2 I Z 2 

B { =Q n + 321.2 - 16.7-^ + 0.218 1 -j- 



(18) 



Note that neither the angular momentum nor the ex- 
citation energy of the nucleus are taken into account 
in the estimate of the fission barriers. 

2) Zj > 89. For heavy fissioning nuclei with 
Zj > 89 GEM follows the RAL model 12 ' 13 and does 
not calculate at all the fission width Tf and does 
not use Eq. (13) to estimate the fission probability 
Pf. Instead, the following semi-empirical expression 
obtained by Atchison 12 ' 13 by approximating the experi- 
mental values of r„/T/ published by Vandenbosch and 
Huizenga 50 is used to calculate the fission probability: 



log(T n /T f ) = C(Zi)(Ai - A (Zi)), 



(19) 



where C(Z) and Ao(Z) are constants dependent on 
the nuclear charge Z only. The values of these 
constants are those used in the current version of 
LAHET 14 and are tabulated in Table 5 (note that 
some adjustments of these values have been done since 
Atchison's papers 12 ' 13 were published). 
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Table 5. C(Z) and A (Z) values used in GEM 



z 


C(Z) 


A (Z) 


89 


0.23000 


219.40 


90 


0.23300 


226.90 


91 


0.12225 


229.75 


92 


0.14727 


234.04 


93 


0.13559 


238.88 


94 


0.15735 


241.34 


95 


0.16597 


243.04 


96 


0.17589 


245.52 


97 


0.18018 


246.84 


98 


0.19568 


250.18 


99 


0.16313 


254.00 


100 


0.17123 


257.80 


101 


0.17123 


261.30 


102 


0.17123 


264.80 


103 


0.17123 


268.30 


104 


0.17123 


271.80 


105 


0.17123 


275.30 


106 


0.17123 


278.80 



In this approach the fission probability Pf is indepen- 
dent of the excitation energy of the fissioning nucleus 
and its angular momentum. 

B.2. Mass Distribution. The selection of the mass 
of the fission fragments depends on whether the fission 
is symmetric or asymmetric. For a pre-fission nucleus 
with Zf/Ai < 35, only symmetric fission is allowed. 
For ZfjAi > 35, both symmetric and asymmetric fis- 
sion are allowed, depending on the excitation energy of 
the fissioning nucleus. No new parameters were deter- 
mined for asymmetric fission in GEM. 

For nuclei with Zf/Ai > 35, whether the fission 
is symmetric or not is determined by the asymmetric 
fission probability P asy 



4870e-°- 36g 
l + 4870e-°- 36E ' 



(20) 



B.2.a. Asymmetric fission. For asymmetric fis- 
sion, the mass of one of the post-fission fragments 
A\ is selected from a Gaussian distribution of mean 
Af = 140 and width <jm = 6.5. The mass of the 
second fragment is A2 = Ai — A\. 
B.2.b. Symmetric fission. For symmetric fission, 
Ai is selected from the Gaussian distribution of mean 
Af = Ai/2 and two options for the width <tm as de- 
scribed below. 

The first option for choosing a M is the original 
Atchison approximation: 



cm 



3.97 + 0.425(£ - B f ) - 0.00212(£ - B f ) 2 , 
25.27, 

(21) 



for (E — Bf) below or above 100 MeV, respectively. 
In this expression all values are in MeV and the fission 
barriers Bf are calculated according to Eq. (18) for 
nuclei with Z t < 88. For nuclei with Z- L > 88, the 
expression by Neuzil and Fairhall 51 is used: 



B f = C-0.36(Zf/Ai), 



(22) 



where C = 18.8,18.1,18.1, and 18.5 [MeV] for odd- 
odd, even-odd, odd-even, and even-even nuclei, respec- 
tively. 

The second option in GEM for ou (used here) was 
found by Furihata 37 ' 38 as: 

°m = C^Zf/Arf + Ci(Zf/Ai) + C 5 (E - B f ) + C . 

(23) 

The constants C 3 = 0.122, C" 4 = -7.77, C 5 = 
3.32 x 10~ 2 , and C 6 = 134.0 were obtained by fit- 
ting with GEM the recent Russian collection of ex- 
perimental fission-fragment mass distributions 52 . In 
this expression, the fission barriers Bf by Myers and 
Swiatecki 53 are used. More details may be found in 
Ref. 38 

B.3. Charge Distribution. The charge distribu- 
tion of fission fragments is assumed to be a Gaussian 
distribution of mean Zf and width oz- Zf is expressed 
as 



Z, = 



Zi + Z[ 



where 



65.5A, 
131 + A? /3 ' 



I = 1 or 2. 



(24) 



(25) 



The original Atchison model uses az = 2.0. An investi- 
gation by Furihata 38 suggests that <rz = 0.75 provides 
a better agreement with data; therefore az = 0.75 is 
used in GEM and in our calculations. 

B.4. Kinetic Energy Distribution. The kinetic 
energy of fission fragments [MeV] is determined by a 
Gaussian distribution with mean e/ and width a ef . 
The original parameters in the Atchison model are: 



e f = O.USZf/A 1 / 3 - 11.4, 



<7 e/ = 0.084e/. 



Furihata's parameters in GEM, which we also use, 



are: 
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f O.UlZf/A 1 / 3 , 
' ' ' 0.104Z?/i4} /3 



24.3, 



(26) 



for Zf/A 1 / 3 < 900 and 900 < Z?/Ai' a < 1800, re- 
spectively, according to Rusanov et al. 52 By fitting the 
experimental data by Itkis et al. 54 , Furihata found the 
following expression for a Cf 



1/3 



CiiZf/A 1 / 3 
C 2 , 



1000) + c 2> 



(27) 



for Zf/A^ 3 above and below 1000, respectively, and 
the values of the fitted constants are C\ = 5.70 x 10~ 4 
and C2 = 86.5 The experimental data used by Furihata 
for fitting are the extrapolated values to the nuclear 
temperature 1.5 MeV by Itkis et al. 54 More details may 
be found in 38 . 

We note that Atchison also has modified his origi- 
nal version using recent data and published 55 an im- 
proved (and more complicated) parameterization for 
many quantities and distributions in his model, but 
these modifications 55 are not yet included either in LA- 
HETorinGEM. 

C. Results from CEM2k+GEM2 
We have merged the GEM2 code with our CEM2k, 
initially keeping all the default options in GEM2. We 
began by concentrating on an analysis of the recent GSI 
measurements in inverse kinematics 6-9 as the richest 
and best data set for testing this kind of model. As 
mentioned above, to understand the role of preequilib- 
rium particle emission, we performed calculations of all 
the reactions we tested both taking into account pree- 
quilibrium particle emission and ignoring it, i.e., going 
directly to GEM2 after the intranuclear cascade stage 
of a reaction described by CEM2k. The size of the 
present paper allows us to present only a few results 
for one reaction measured at GSI, which we choose to 
be p(800 MeV) + Au 9 . Results for other reactions may 
be found in 56 . 

If we merge GEM2 with CEM2k without any mod- 
ifications, the new code does not describe correctly 
the fission cross section (and the yields of fission frag- 
ments) whether we take into account preequilibrium 
emission (see the dashed red line on Fig. 5) or not 
(see the dashed blue line on Fig. 5). Such results were 
anticipated, as Atchison fitted the parameters of his 
RAL fission model when it was coupled with the Bertini 
INC 57 which differs from our INC. In addition, he did 
not model preequilibrium emission. Therefore, the dis- 
tributions of fissioning nuclei in A, Z, and excitation 
energy E* simulated by Atchison differ significantly of 
the distributions we get; as a consequence, all the fis- 
sion characteristics are also different. 
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Figure 5. Comparison of the experimental 9 mass and 
charge distributions of the nuclides produced in the 
reaction p(800 MeV) + Au (circles) with different 
calculations. The black solid lines show our calcula- 
tions with CEM2k without contributions from fission. 
The dashed red and blue lines show results found by 
merging CEM2k with GEM2 without any modifications 
when preequilibrium emission is (red lines) or is not 
(blue lines) included. Solid color lines show results from 
CEM2k+GEM2 with a modified a/: red lines are for 

,CBM I „RAL _ 



the case with preequilibrium emission ( 



7 



1.0848) and blues lines show the results without pree- 



quilibrium emission (c 



CEM I RAL 
f l a f 



0.97). 
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Furihata used GEM2 coupled either with the Bertini 
INC 57 or with the ISABEL 58 INC code, which also dif- 
fers from our INC, and did not include preequilibrium 
particle emission. Therefore the real fissioning nuclei 
simulated by Furihata differ from the ones in our sim- 
ulations, and the parameters adjusted by Furihata to 
work the best with her INC should not be the best 
for us. To get a good description of the fission cross 
section (and fission-fragment yields) we need to mod- 
ify at least one parameter in GEM2, namely to adjust 
the level density parameter a/ to get the correct fis- 
sion cross section (see Eq. (17)), in the case of fis- 
sioning nuclei with Z < 88 (pre-actinides), and the 
parameter C(Z) (see Eq. (19) and Tab. 5) for fission- 
ing nuclei with Z > 88 (actinides). From the dashed 
lines on Fig. 5 we see that we need to enlarge a/ in 
our code to get a proper fission cross section when 
we include preequilibrium emission (the excitation en- 
ergy of our fissioning nuclei and their A and Z are 
smaller than provided by the Bertini or ISABEL INC 
without preequilibrium), and we need to decrease a/ 
in the case without preequilibrium. By increasing a/ 
by 1.0848 compared with the original RAL and GEM2 
value (a < f EM /af AL = 1.0848), we are able to repro- 
duce correctly the fission cross section when we take 
into account preequilibrium emission (below, we label 
such results as "with Prec"). In the case with no pree- 
quilibrium emission, a proper fission cross section is 
obtained for a c f EM /af AL = 0.97 (we label such re- 
sults as "no Prec"). We choose these values for a/ for 
all our further calculations of this reaction and do not 
change any other parameters. 

The solid lines in Fig. 5 show results with these val- 
ues of a,f. One can see that the "no Prec" version pro- 
vide a good description of both the mass and charge 
distributions and agrees better with the data for these 
characteristics than the "with Prec" version (that is not 
true for isotopic distributions of individual elements, as 
we show below). The "with Prec" version reproduces 
correctly the position of the maximum in both A and Z 
distributions and the yields of fission fragments not too 
far from these maximums, but the calculated distribu- 
tions are narrower than the experimental ones. This is 
again because both Atchison and Furihata fitted their 
A and Z distributions using models without preequi- 
librium emission, which provide higher values for the 
excitation energy, A, and Z of fissioning nuclei. This 
means that to get a good description of A and Z dis- 
tributions for fission fragments using GEM2 in CEM2k 
"with Prec", we would need to modify the A and Z 
distributions of fission fragments in GEM2 (see Sec. 
B.2 and B.3), making them wider. This would take us 
beyond the scope of the present work and here we do 
not vary any more parameters than we have already 
discussed. 

Fig. 6 shows the GSI measurements 9 of the A and Z 
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Figure 6. Comparison of the experimental 9 mass and 
charge distributions of the kinetic energy of the nu- 
clides produced in the reaction p(800 MeV) + Au 
(circles) with out CEM2k+GEM2 calculations: "with 
Prec" results are shown by red lines, "no Prec" results 
are shown by blue lines. 

distributions of the kinetic energy of products from the 
same reaction compared with our CEM2k+GEM2 cal- 
culations both with and without preequilibrium emis- 
sion. Both versions of our calculations are in reasonable 
agreement with the data, with slightly better agree- 
ment for the version "no Prec". These results suggest 
to us that a small adjustment of kinetic-energy dis- 
tribution parameters in GEM (see Sec. B.4) may be 
required for a better description of the data with our 
CEM2k+GEM2 code. 

Mass and charge distributions of the yields or kinetic 
energies of the nuclides produced show only general 
trends and are not sensitive enough to the details of 
a reaction. It is much more informative to study the 
characteristics of individual nuclides and particles pro- 
duced in a reaction. Fig. 7 shows a comparison of the 
experimental data on production yields of twelve sep- 
arate isotopes with Z lying from 20 to 80 from the 
same reaction measured at GSI 9 with our calculations 
using both the "with Prec" (upper plot) and "no Prec" 
(lower plot) versions. 
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Figure 7. Experimental 9 mass distributions of the cross 
sections of twelve isotopes with the charge Z from 20 
to 80 compared with our CEM2k+GEM2 calculations. 
"With Prec" results are shown on the upper plot, while 
"no Prec" results are shown in the lower one. 

The agreement (or disagreement) of our calculations 
with these data is different from what we have for the 
integral A or Z distributions in Figs. 5 and 6: We see 
that for the isotopes produced in the spallation region 
(not too far from the target) and for fission fragments 
in the region with the maximum yield, the version "with 
Prec" agree much better with the data than the ver- 
sion "no Prec" . Only for production of isotopes at the 
border between spallation and fission and between fis- 
sion and fragmentation does the version "with Prec" 
underestimates the data, due to too narrow A and Z 
distributions in the simulation of fission fragments, as 
we discussed previously. The "no Prec" version agrees 
better with the data in these transition regions but are 
in worse agreement for isotopes both in the spallation 
region and in the middle of the fission region. We con- 
clude that if a model agrees well with some A or Z 
distributions it does not necessarily mean that it also 
describes well production of separate isotopes. In other 
words, integral A and Z distributions are not sensitive 
enough to develop and test such models, a practice 
which has often used in the literature. 

It is more difficult for any model to describe cor- 
rectly the energy dependence for the production cross 
sections of different isotopes, i.e.. excitation functions. 
We calculated using both the "with Prec" and "no 
Prec" versions of CEM2k+GEM2 all the excitation 
functions for the same reaction, p + Au, for proton en- 
ergies from 10 MeV to 3 GeV and compared our results 
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Figure 8. Excitation functions for the production of 

196 Au _ 193 Hgi 194 Hgj and 195 Hg from p+ 197 Au Re _ 

suits by CEM2k+GEM2 "with Prec" are shown by solid 
lines and "no Prec" by dashed lines. Experimental data 
(symbols) are from our LANL compilation (T16 Lib) 59 
and are available from the authors upon request. 

with all available data from our compilation referred to 
here as T-16 Library ("T16 Lib") 59 . Only several typ- 
ical examples from our comparison are shown below. 
Figs. 8 and 9 show two examples of excitation functions 
for the production of several isotopes in the spallation 
region. One can see a not too good but still reasonable 
agreement of both calculations with the data. Never- 
theless, the version "with Prec" reproduces these ex- 
perimental excitation functions better that the version 
"no Prec". Similar results were obtained for excita- 
tion functions of many other isotopes in the spallation 
region. 

Figs. 10 and 11 show two examples of excitation 
functions for the production of fission fragments. We 
see that merging CEM2k with GEM2 allows us to rea- 
sonably describe yields of fission fragments, while in the 
old standard CEM2k we do not have any fission frag- 
ments and are not able to describe such reactions at all. 
We see that as shown in Figs. 5 and 7 for a single pro- 
ton energy of 800 MeV, the "with Prec" version agrees 
better with the data in the whole energy region for the 
production of most of the fission fragments. Only on 
the border between fission and fragmentation regions 
( 54 Mn and 60 Co in Fig. 10) does the "no Prec" version 
agree much better with the data than the "with Prec" 
version; the reason for this we have already discussed. 

In Fig. 12 we show examples of excitation functions 
for the production of light fragments, in the fragmenta- 
tion region, that are produced in CEM2k+GEM2 only 
via evaporation (the contribution to the yield of these 
isotopes from fission or deep spallation is negligible). 
We see that with the "no Prec" version, GEM2 repro- 
duces correctly the yields of light fragments 7 Be and 
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Figure 10. The same as Fig. 8 but for the production 
of 54 Mn, 60 Co, 74 As, 82 Br, 84 Rb, and 88 Y. Note that 
while some minor contribution to the production of 
these isotopes from deep spallation processes may be 
present, all of them are produced mainly via fission. 



Figure 9. The same as Fig. 8 but for the produc- 
tion of !53Gd cum , 168 Tm, I89| r c " m , 190 lr, 192 lr, and 
194 Au. Note that for 153 Gd, only the cumulative yield 
is shown, where measured and calculated cross sections 
contain contributions not only from a direct produc- 
tion of 153 Gd ("independent yield"), but also from all 
its decay-chain precursors. For 189 lr, both cumulative 
and independent cross sections are shown. 
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Figure 11. The same as in Fig. 8 but for the production 
of 95 Nb, 96 Tc, 102 Rh, and 110 Ag. 
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Figure 12. The same as Fig. 8 but for the produc- 
tion of 7 Be, 18 F CMm , 24 Ne cmm , 22 Na cum , 24 Na CMm , 
and 28 Mg CMm . These fragments are produced in 
CEM2k+GEM2 only via evaporation. 

18 F, and not so well the excitation functions for heavier 
fragments 22 Na and 24 Ne. With increasing mass of the 
fragment, the calculations progressively underestimate 
their yields (e.g., for 28 Mg, the calculated excitations 
function is more than an order of magnitude below 
the data). The version "with Prec" strongly underes- 
timates the yields of all these fragments, and this is 
again not surprising, as Furihata developed her model 
and fitted all parameters without taking into account 
preequilibrium processes. Undeniably, the parameters 
determining the yields of evaporated fragments in GEM 
(inverse cross sections and Coulomb barriers) could be 
adjusted to get a good agreement with the data for 
the yields of light fragments with the version "with 
Prec". This is not an aim of our present work and we 
will not do this here. Even if we were to do this, we 
expect in advance to get similar results as we got for 
the "no Prec" version: It would be possible to describe 
correctly the yields of light fragments but not of heavy 
fragments like 24 Na and 28 Mg. To describe such heavy 
fragments the model would need to be improved fur- 
ther, by considering other mechanisms for heavy frag- 
ment production in addition to the evaporation process 
taken into account by GEM2. 
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Figure 13. The same as Fig. 8 but for the production 
of n, p, d, t, 3 He, and 4 He. The complex particles are 
produced in CEM2k+GEM2 via evaporation and pree- 
quilibrium emission; n and p are also produced during 
the cascade stage. 
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Finally, Fig. 13 shows the excitation functions for 
emission of nucleons and complex particles up to a for 
this reaction. Note that the data for these excitation 
functions are not so extensive and precise as we have 
for heavier products: many data points were obtained 
by integration (plus extrapolation) of the spectra of 
particles measured only at several angles and only for 
a limited range of energy. But even from a comparison 
with these sparce and imprecise data we see that the 
"with Prec" version describes these excitation function 
much better than the "no Prec" version. This is an 
expected result as the high Coulomb barriers for heavy 
nuclear targets oppose evaporation of low energy com- 
plex particles and the main contribution to their yields 
comes from preequilibrium emission from highly excited 
pre-compound nuclei. 

Besides yields of reaction products, reliable models 
should also describe their spectra. Data on particle 
spectra are important for shielding calculations and 
other applications, in addition to the scientifc interest 
in understanding nuclear reactions. For the reaction 
discussed here, p + Au, there are measurements by 
Bertrand and Peelle for proton and complex-particle 
spectra at incident proton energies of 29, 39, and 62 
MeV 60 . We compare calculated angle-integrated en- 
ergy spectra for p, d, t, 3 He, and 4 He with the data at 
62 MeV in Fig. 14. 

We see that the "with Prec" version (black solid his- 
tograms) describes a good part of the complex particle 
emission spectra, though the tails of calculated spectra 
are significantly below the data (we discuss in the next 
section an attempt to improve the description of com- 
plex particle spectra in CEM2k). The "no Prec" ver- 
sion (dashed red and green histograms) fails completely 
to describe these spectra, as it contains contribution to 
complex particle emission only from evaporation, and 
evaporation spectra do not extend significantly above 
particle energies of 20 MeV. Even in the low-energy 
evaporation region, the evaporation component for the 
3 He spectra (shown by a green dashed histogram in 
Fig. 14) is more than an order of magnitude below the 
data 60 . Let us mention that this is a problem not only 
of GEM2 but also for all other similar models where 
preequilibrium emission is not modeled. What is more, 
some of the evaporation models used in the literature, 
like the GSI evaporation model by Schmidt et al. 61 , 
which is used in conjunction with the Liege INC by 
Cugnon et al. 15 , evaporate only n, p, and 4 He, and 
do not include evaporation of d, t, and 3 He providing 
for them no yield. 

For completeness sake, we show here also an ex- 
ample of results from a calculation with the merged 
CEM2k+GEM2 code of a reaction on an actinide, 
p(190 MeV) + 232 Th. This reaction was recently mea- 
sured by Duijvestijn et al. 62,63 , and at a similar energy 
of 200 MeV, by Titarenko et al. 64 To get for actinides 
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Figure 14. Angle-integrated energy spectra of p, d, t, 
3 He, and 4 He from 61.5 MeV protons on 197 Au. Cal- 
culations from the merged CEM2k+GEM2 code "with 
Prec" are shown by solid black histograms and "no 
Prec", by dashed red (and green, for 3 He) histograms. 
Experimental data (symbols) are from Bertrand and 
Peelle 60 . 
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a proper fission cross section, we need to adjust in 
GEM2 the parameters C(Z) (or, also A (Z)) in Eq. 
(19), as they were fitted by Atchison to work the best 
with Bertini's INC and we have in CEM2k our own 
INC. As mentioned above, for actinides, Eq. (13) is 
not used in GEM2 and a/ is not used in any calcula- 
tions, therefore we do not need to adjust ct//a n , for 
fissioning nuclei with Z > 88. We found that to get 
with CEM2k+GEM2 a fission cross section in agree- 
ment with the data for our reaction we need to use 
C{Z) CEM /C{Z) RAL = 8.50 for the "with Prec" ver- 
sion and C(Z) CEM /C(Z) RAL = 1.77 for "no Prec" 
(while we use a c f EM /af AL = 1.0848 which we fit- 
ted for the reaction p(800 MeV) + Au, "with Prec"). 
This is the only parameter we fitted for this reaction. 
Nevertheless, we should mention that for reactions 
on actinides at intermediate or high energies, the pa- 

CEM j„RAL should a | s0 De fj ttec j a | ong w j tn 



rameter a 



7«f 



C(Z) CEM /C{Z) RAL . In some simulated events sev- 
eral protons can be emitted at the cascade and preequi- 
librium stages of the reaction, as well as at the evap- 
oration stage, before the compound nucleus actually 
fissions, and the charge of the fissioning nucleus can 
have Z < 88, even when the initial charge of the tar- 
get has Z > 88. At the same time, for Z < 88, due to 
charge exchange reactions, the charge of the fission- 
ing nucleus may exceed 88, so that we would need to 
fit as well C(Z) CEM /C(Z) RAL . This is a peculiar- 
ity of treating the fission probability Pf differently for 
the elements above and below Z = 89 in the Atchison 
model. 

Fig. 15 shows mass distributions of products from 
p(190 MeV) + 232 Th calculated with both versions of 
CEM2k+GEM2 compared to the available experimen- 
tal data 62 ' 64 . 

We need to mention that these data are not as good 
for testing and developing models as are the GSI data 
measured in inverse kinematics for the p + Au reaction 
discussed above: All the data shown in Fig. 15 were 
obtained by the 7-spectrometry method. Only some 
of the produced isotopes were measured, and most of 
the data were measured for the cumulative yields. To 
get the "experimental" A-distribution, we summed for 
each A the available data taking care to not sum the in- 
dividual cross sections already included in some cumu- 
lative yields; but the resulting A-distribution is still not 
complete, as many isotopes were not measured. This 
means that some theoretical values can be above the 
experimental data (where some isotopes were not mea- 
sured) without necessarily implying disagreement be- 
tween calculations and measurements. What is more, 
many of Duijvestijn's data tabulated in 62 differ sig- 
nificantly from the same data tabulated in an earlier 
publication 63 for some isotopes: e.g., for 72 Zn the dif- 
ference is a factor of 8.03 (!), possibly due to a misprint 
in 63 . We chose for our comparison Duijvestijn's data 
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Figure 15. Mass distribution of nuclides from p (190 
MeV) + 232 Th. The 190 MeV data are from 62 ; the 200 
MeV data are by Titarenko et al. , data for several Ac 
isotopes from p (150 MeV) + 232 Th are from the T-16 
compilation 59 . Our "with Prec" results are shown by 
the solid blue line and the "no Prec" results are shown 
by the dashed red line. 

tabulated in 62 as more reliable than the earlier tabu- 
lation published in 63 . (We thank Dr. Duijvestijn for 
clarifying this point.) 

One can see that after the parameter C(Z) is ad- 
justed, both the "with Prec" and "no Prec" versions 
of CEM2k+GEM2 describe equally well the mass dis- 
tribution of the products from this reaction; therefore 
it is not possible to choose between either of the ver- 
sions from a comparison with these not very informa- 
tive data. It is more useful to compare calculations 
with the yields of individually measured isotopes. In 
Fig. 16, we compare our calculations with all cross sec- 
tions measured by Duijvestijn 62 for each nuclide sep- 
arately, where we can compare our results with the 
data (we do not include in our comparison the nu- 
clides measured only either in their isomer or ground 
states, as our model does not provide such information: 
CEM2k+GEM2 provides only yields for the sum of iso- 
tope production cross sections both in their ground and 
excited states). We see that on the whole, the "with 
Prec" version agrees better with most of the individu- 
ally measured cross sections than the "no Prec" version 
and for many of the measured isotopes the disagree- 
ment is less than a factor of two. Nevertheless, for sev- 
eral isotopes like 72 Ga, 96 Tc, and 1 I, we see some big 
disagreements. For comparison, we also show in Fig. 
16 calculations by the phenomenological code YIELDX 
of Silberberg, Tsao, and Barghouty 23 and with the 
phenomenological code CYF by Wahl 65 often used in 
applications. We see that both these phenomenologi- 
cal systematics completely fail to describe the produc- 
tion of all isotopes from this reaction, indicating that 
we cannot rely on phenomenological systematics and 
must develop reliable models to be used in applications. 
We note that the agreement of our CEM2k+GEM2 re- 
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suits with the measured product yields is different for 
different reactions. So, in Ref. 56 , for the reaction p 
(100 MeV) + 238 U, we got an almost perfect agree- 
ment of calculations with the data: for this particular 
reaction, both "with Prec" and "no Prec" results by 
CEM2k+GEM2 almost coincide with the available data 

(see details in Ref. 56 ). 

To summarize this section, merging CEM2k with 

GEM2 allows us to describe reasonably well many fis- 
sion and fragmentation reactions in addition to the 
spallation reactions already described well by CEM2k. 
Some reactions, like production of fission fragments 
at the borders between fission and fragmentation or 
between fission and emission of heavy fragments like 
Na and Mg are poorly described by CEM2k+GEM2 in 
its current version. This disagreement does not dis- 
courage us; the results of the present work suggest 
that some of the fission and evaporation parameters of 
GEM2 can be adjusted to get a much better description 
of all reactions. This lends credibility to such an ap- 
proach. There is one more drawback of this approach 
to be mentioned: considering evaporation of up to 66 
particles in GEM becomes extremely time consuming 
when calculating reactions with heavy targets at high 
incident energies. But even this disadvantage may be 
mitigated by the performance of modern computers. 
We have nevertheless some more serious doubts about 
the current version of GEM2 related to its lack of self- 
consistency, e.g.: 

1) using different, not physically related parameteri- 
zations for inverse cross sections and Coulomb barriers 
for different particles and fragments; 

2) using different level density parameters for the 
same compound nuclei when calculating evaporation 
and estimating fission probability from the widths of 
neutron evaporation and fission; 

3) different, and purely phenomenological treat- 
ments of fission for pre-actinide and actinide nuclei; 

4) not taking into account at all the angular mo- 
mentum of compound and fissioning nuclei; 

5) rough estimations for the fission barriers and level 
density parameters, etc. 

This means that an approach like GEM2 can in prin- 
ciple be used to describe fission and evaporation of 
particles and fragments heavier that He after the INC 
and preequilibrium parts of CEM2k and other models, 
but it should be considerably improved striving first to 
progressively incorporate better physics, and only after 
that looking on agreement with the data. 

The results of the present work and from 56 show 
that on the whole, the merged CEM2k+GEM2 code 
agrees better with most of tested experimental data 
when we take into account the preequilibrium emis- 
sion of particles, than when we neglect completely 
preequilibrium processes. But there is still an open 
question to be solved here; we have had some indica- 
tions for many years that CEM accounts for too many 
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Figure 16. Detailed comparison 

between experimental 62 (black filled squares) and cal- 
culated cross sections of individual (upper plot) and cu- 
mulative (middle and bottom plots) reaction products. 
Our CEM2k+GEM2 "with Prec" results are shown by 
connected blue squares and the "no Prec" results are 
shown by connected red circles. For comparison, pre- 
dictions by the phenomenological systematics YIELDX 
of Silberberg, Tsao, and Barghouty 23 are shown by 
connected indigo diamonds, and results calculated by 
Dr. W. B. Wilson with the phenomenological code CYF 
by Wahl 65 are shown with connected green diamonds. 
We thank Dr. Wilson for providing us with results of his 
calculations with the CYF code included in this figure. 
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preequilibrium particles, at least at energies above the 
pion-production threshold. To solve this problem, as 
a "zero-step" approximation, in our original CEM2k 
version 10 we neglected in the exciton model the tran- 
sitions that decrease or do not change the number of 
excitons A„ = —2 and A n = 0, shortening in this way 
the preequilibrium stage of reactions. We do not like 
this approach as it is arbitrary, "ad hoc", even though 
this "never come back" approximation is used in some 
popular codes 14,66 . In the present work we removed 
this arbitrary condition in CEM2k and the "with Prec" 
version takes into account all the preequilibrium tran- 
sitions A n = +2, 0, and -2, making the preequilibrium 
stage of a reaction longer and increasing the number 
of emitted preequilibrium particles. The results of the 
present work indicate to us once again that we need to 
take into account the preequilibrium stage in reactions, 
but we need less particle emission than we currently 
calculate at this stage. In 56 , we explore an alternative 
way to shorten the preequilibrium particle emission in 
CEM, based on Ref. 67 , but this work is incomplete so 
we do not include these results here. 

IV. EXTENSION OF CEM2K 
In this section we present an alternative way of de- 
scribing fission-fragment, complex-particle and light- 
fragment emission with CEM2k, by extending and de- 
veloping further our model itself, without using inde- 
pendent evaporation models developed by other au- 
thors. 

One of the unsolved problems in all versions of the 
CEM code is using the Dostrovsky et. al. 39 approxi- 
mations for the inverse cross sections (Eq. (2) with the 
parameters shown in Tab. 3) both for preequilibrium 
and evaporation. As already mentioned, the Dostro- 
vsky et. al. 39 formula is simple and allows one to cal- 
culate the integral in Eq. (1) analytically. This is the 
main reason we have kept Dostrovsky's formula in all 
our previous versions of the CEM code, even though 
we know it is not reliable enough. We develop and 
incorporate into CEM2k our own approximation for in- 
verse cross sections. As a first step, we collected ex- 
perimental data on absorption or fusion cross sections 
for all target nuclei and for all particles and fragments 
for which we were able to find data, i.e., just the "in- 
verse" cross sections when we deal with particle emis- 
sion. Then, we collected from the literature theoretical 
estimations and systematics for the inverse cross sec- 
tions and compared these systematics with the data to 
find which systematics better describe the data. An 
example from this study is shown in Fig. 17, where we 
compare available data on inverse cross sections for n, 
p, d, t, 3 He, and 4 He on 27 Al with the approxima- 
tion from Ref. 39 , NASA systematics by Tripathi, Cu- 
cinota, and Wilson 33 , a parameterization by Kalbach 34 , 
an approximation by Tang, Srinivasan, and Azziz 69 , 



and the results calculated with the phenomenological 
code CROSEC by Barashenkov and Polanski 70 . Simi- 
lar figures have been plotted for other target nuclei for 
which we were able to find experimental data. One can 
see a big disagreement between Dostrovsky et al.'s 39 
approximation and the data. 

Probably, the reason why Dostrovsky et. al.'s 39 for- 
mula works quite well in many evaporation models to 
describe emission of neutrons and protons from many 
reactions is because it overestimates both and in simi- 
lar ways the neutron and proton inverse cross sections 
(correspondingly, their width T n and T p ) at energies 
above about 20 MeV. The emission of n and p are 
the main channels of most evaporation processes, and 
emission of particles are simulated in evaporation codes 
using only the ratios r,-/£) . Tj but not the absolute 
values of particle widths, Tj. Dostrovsky formula un- 
derestimates significantly the inverse cross sections for 
the complex particles, and this is one of the reasons 
why all evaporation models that use Dostrovsky's for- 
mula underestimate the yields of complex particles. An 
additional problem introduced by the overestimation of 
the nucleon emission probabilities occurs when fission 
is considered. A realistic value of the fission barrier 
and the saddle-point level density will give too small 
a fission probability when compared to a too-large nu- 
cleon evaporation probability. This then leads to the 
empirical determination of unphysically large values of 
the level density parameter, etc. 

From Fig. 17 and similar figures for other nuclei 
which we have made, we determined that the over- 
all best agreement with experimental data is achieved 
with the NASA systematics by Tripathi, Cucinota, and 
Wilson 33 ; therefore we chose it as the basis to calculate 
inverse cross sections in CEM2k. Nevertheless, this sys- 
tematics does not reproduce correctly the inverse cross 
sections for neutrons at energies below about 10 MeV. 
We address this problem the following way: We cal- 
culate with the NASA systematics 33 the inverse cross 
sections for all charged particles and for neutrons with 
energies above T" Lax , where a maximum in the neu- 
tron inverse cross section is predicted by 33 , and with 
the systematics by Kalbach 34 , for neutrons with ener- 
gies below T™ aa; . We tabulated the values of T™ ax for 
each nucleus and renormalize the neutron inverse cross 
sections calculated with the systematics by Kalbach 34 
so that they coincide at T™ aa: with the NASA values 33 , 
providing a continuity of the neutron inverse cross sec- 
tions at this energy. 

A detailed description of the NASA and Kalbach sys- 
tematics may be found in 33,34 and references therein. 
For completeness sake, we outline here only their basic 
ideas and formulas. 

The NASA systematics 33 is a universal parameter- 
ization for any systems of colliding nuclei, therefore 
can be used to calculate inverse cross sections not only 
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for nucleons and complex particles but also for heavier 
fragments. It uses the following form for the reaction 
cross sections (inverse cross sections, in our case): 

o R = nr 2 (A p /3 + A T /3 + S E f(l - V c /E cm ), (28) 

where Ap and At are the projectile and target mass 
numbers, respectively, r = 1.1 fm is energy inde- 
pendent, E cm is the total center-of-mass kinetic en- 
ergy in MeV, Se is an energy dependent parameter de- 
scribed below, and Vc [MeV] is the energy-dependent 
Coulomb barrier defined as: 



V c = 1.4AZ P Z T /R. 



(29) 



Here, 



R = r P + r T + 1.2(4 /3 + A 1 ^) / E 1 ^ , (30) 
with (i = P, T) 



n = 1.29(r;) T 



(31) 



with the root-mean-square radius, (rj) rms , obtained 
directly from experimental data. 

There is an energy dependence in the reaction cross 
section at intermediate and higher energies mainly due 
to two effects — transparency and Pauli blocking. This 
is taken into account in Se which is given by 



S E = 1.855 + 0.16- 



S 



E, 



1/3 



+ 0.91 



{A t - 2Z T )Zj 



At A 



(32) 



where S is the mass asymmetry term and is given by 

(33) 



5 = 4 /3 4 /3 /(4 /3 +4 /s ) 



and is related to the volume overlap of the collision sys- 
tem. The last term on the right hand side of Eq. (32) 
accounts for the isotopic dependence of the reaction 
cross section. 

The term Ce is related to the transparency and Pauli 
blocking and is given by 



C E = D (1 - exp(-£/40)) 
x cos(0.229£ - 453 ). 



0.292 exp(-£/792) 
(34) 



Here D is related to the density dependence of the 
colliding system scaled with respect to the density of 
the C+C system, i.e.: 



D = 1.75(pA P + PA T )/{pA c + PAc)- 



(35) 



The density of a nucleus is calculated in the sharp- 
surfaced sphere model and for a nucleus of mass Ai is 
given by 



PAi = Ai/-TTrf, 



(36) 



where the radius of the nucleus is defined in Eq. (31). 
The physics related to the constant D is to simulate 
the modifications of the reaction cross sections due to 
Pauli blocking. This effect helps present a universal 
picture of the reaction cross section 33 . 

At lower energies (below several tens of MeV) where 
the overlap of interacting nuclei is small (and where the 
Coulomb interaction modifies the reaction cross section 
significantly) the modifications of the cross sections 
due to Pauli blocking are small, and gradually play 
an increasing role as the energy increases, since this 
leads to higher densities where Pauli blocking gets in- 
creasingly important. For proton-nucleus interactions, 
where there is not much compression effect, a single 
constant value of D = 2.05 gives good results for all 
proton-nucleus collisions. For alpha-nucleus collisions, 
where there is a little compression, the best value of D 
is given by 33 

D = 2.77 - 8.0 x 10~ 3 A T + 1.8 x 1Q- 5 A T 

- 0.8/(1 + exp(250- £)/75). (37) 

For lithium nuclei because of the " halos" , compression 
is less and hence the Pauli blocking effect is less impor- 
tant and a reduced value of D/3 gives better results 
for the reaction cross sections at the intermediate and 
higher energies. 

Note that for proton-nucleus collisions this method 
of calculating the Coulomb energy underestimates its 
value for the very light closed shell nuclei of alpha and 
carbon, and these should be increased by a factor of 
27 and 3.5 respectively for a better fit 33 . 

Kalbach's approximation for inverse cross sections 34 
is based on optical-model reaction cross sections using 
an empirical parameterization of Chatterjee, Murthy, 
and Gupta 73 fitted for different particles using differ- 
ent optical potentials, with some additional modifica- 
tions and described in 34 . For instance, the reaction 
cross section ap (inverse cross section, in our case) for 
neutrons (used here) is written as 



v/e, 



(38) 



where A, fi, v are mass-dependent parameters and e is 
the neutron laboratory energy in MeV. As mentioned 
in 73 , the addition of the linear term Ae in Eq. (38) 
greatly improves the fit as compared with the Dostro- 
vsky et al. approximation 39 . Separating out the energy 
dependence as in Eq. (38), the dependence of A, /x, v 
on target mass number A, was obtained empirically 



as 



73. 



A = AoA- 1/3 + A l7 

/i = ^ 1/3 + Mi^ 2/3 , 
^ = ^ 4 / 3 + ^A 2 / 3 + z, 2 . 



(39) 
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Figure 17. Comparison of experimental data on absorption cross sections for n, p, d, t, 3 He, and 4 He on 27 AI 
with the Dostrovsky et. al. 39 approximation, the NASA systematics by Tripathi, Cucinota, and Wilson 33 , a 
parameterization by Kalbach 34 , an approximation by Tang, Srinivasan, and Azziz 69 , and results calculated with 
the phenomenological code CROSEC by Barashenkov and Polanski 70 . Most of the data points (circles, either 
without a special label or labeled as "BAR93") shown here are from the compilation by Barashenkov 35 ; several 
data points are by Auce et al. 71 (triangles labeled as 'AUC94"), and a data point for 4 He is measured by Dubar 
et al. 72 (the square labeled as "DUB89" ). 
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The seven parameters A , Ai, Mo, Mi, v§, v\, and 
vi in Eq. (39) were fitted and tabulated by Chatterjee, 
Murthy, and Gupta 73 for different optical potentials 
(for different incident particles from n up to a), with a 
priority to the potential by Mani et al. 74 for neutrons. 
In 1982, the original 73 neutron parameterization was 
modified by S. K. Gupta to include an effective bar- 
rier of 2.4 MeV (see details in 34 ), and is now used in 
the Kalbach approximation for neutron inverse cross 
sections 34 with the following values (different from 
the ones of Ref. 73 ) for the parameters in Eq. (38): 
A = 12.10, Ai = -11.27, mo = 234.1, mi = 38.26, 
vq = 1.55, v x = -106.1, and v 2 = 1280.8. In addi- 
tion, Kalbach has changed somehow arbitrarily 34 the 
neutron "barrier" from 2.4 MeV to 0.5 MeV to raise 
the neutron emission cross section for emission ener- 
gies below about 1.5 MeV. Also, after the choice of 
the neutron potential of Mani et al. 74 was made in 
1963, new data became available, and the earlier com- 
parisons with measured non-elastic cross sections have 
been extended to include additional targets. This led 
to a normalization factor for the parameterized optical 
model cross section having the form 34 : 

R m _ J 0.7 + 0.3,4/40 for^<40, 
Rn{A) -\l.O for A > 40. (40) 

Thus, the current version of the Kalbach approximation 
for neutron inverse cross sections used here reads like 
&inv = &r x Rn{A) with ctr and R n {A) defined by 
Eqs. (38-40) and the seven parameters listed above. 

An example of inverse cross sections calculated with 
our method compared to experimental data for n, p, d, 
t, 3 He, and 4 He on 12 C, 40 Ca, 56 Fe, and 238 U is shown 
in Fig. 18. 

For comparison, calculations using the NASA sys- 
tematics by Tripathi, Cucinota, and Wilson 33 , the 
Kalbach parameterization 34 , and the Dostrovsky et. 
al. 39 approximation are shown as well. One can see 
that our hybrid approach combining the NASA sys- 
tematics by Tripathi, Cucinota, and Wilson 33 and the 
Kalbach parameterization 34 reproduces well the mea- 
sured inverse cross sections for n, p, d, t, 3 He, and 
4 He and agrees much better with the data than the 
approximation by Dostrovsky et. al. 39 used earlier in 
all versions of our CEM code. Similar results are ob- 
tained for other nuclear targets for which we found 
data. 

We incorporated this hybrid approximation for ai nv 
in our CEM2k code and use it to calculate inverse cross 
sections for all particles to be evaporated from a com- 
pound nucleus or emitted at the preequilibrium stage 
of a reaction. The widths for particle emission, Tj, 
are calculated by integrating Eq. (1) (and a similar 
equation for preequilibrium particles; see details in 68 ) 
numerically. 

To be able to describe production of light fragments 



heavier than 4 He, as a first step, we extended the evap- 
oration process including emission of up to 66 different 
particles, the same considered by Furihata in GEM2 
and listed in Tab. 1. Following Furihata, we consider 
evaporation of fragments both in the ground and ex- 
cited states, and calculate the widths for the emission 
of excited fragments and their energy in the same man- 
ner as realized in GEM2 by Furihata and described 
above in Sec. III. A. 

We enlarged the number of fragments that may be 
emitted at the preequilibrium stage, considering the 
possibility of emission of up to 29 different preequilib- 
rium particles and fragments (the first 29 from Tab. 1, 
those having Z < 6). 

As we use now in CEM2k the Tripathi, Cucinota, and 
Wilson systematics 33 for charged-particle inverse cross 
sections that employ the energy dependent approxima- 
tion for the Coulomb barriers defined by Eqs. (29-31), 
we replaced the old routine that calculates Coulomb 
barriers according to Dostrovsky et. al. 39 used in all 
previous versions of CEM code with a new routine that 
calculates Coulomb barriers using Eqs. (29-31). We 
did this for self-consistency of the model: we use the 
same approximation for Coulomb barriers while calcu- 
lating inverse cross sections and in all other parts of 
the code where Coulomb barriers are used. 

In this extended version of the CEM2k, we consider 
also the possibility of "creating" high-energy d, t, 3 He, 
and 4 He by final state interactions among emitted cas- 
cade nucleons outside of the target, using the coales- 
cence model as implemented by Gudima et al. 75 In con- 
trast to most other coalescence models for heavy-ion 
induced reactions, where complex-particle spectra are 
estimated simply by convolving the measured or calcu- 
lated inclusive spectra of nucleons with corresponding 
fitted coefficients (see, e.g. 76 , and references therein), 
CEM2k uses in its simulation of particle coalescence 
real information about all emitted cascade nucleons 
and does not use integrated spectra. CEM2k assumes 
that all the cascade nucleons having differences in their 
momenta smaller than p c and a correct isotopic content 
form an appropriate composite particle. This means 
that the formation probability for, e.g., a deuteron is 

W d (p, t) = j j dp p dp n p c (p p , t)p c (p n , t) 

xS(p p +p n -p)Q( Pc - \p c - m - - n c - m -\), (41) 

where the particle density in momentum space is re- 
lated to the one-particle distribution function / by 

p C (p,t) = J dff c (f,p,t). (42) 

Here, the superscript C shows that only cascade nucle- 
ons are taken into account for the coalescence process. 
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Figure 18. Inverse cross sections for n, p, d, t, 3 He, and 4 He on 12 C, 40 Ca, 56 Fe, and 238 U calculated with our 
routine HYBRID (black solid lines) compared with available experimental data (symbols) from the compilation 
by Barashenkov 35 (circles, either without a special label or labeled as "BAR93"), by Auce et al. 71 (labeled 
as "AUC94"), and by Dubar et al. 72 (labeled as "DUB89"). For comparison, calculations using the NASA 
systematics by Tripathi, Cucinota, and Wilson 33 are shown by magenta dot-dashed lines, results from the Kalbach 
parameterization 34 are shown with red dashed lines, and predictions by the Dostrovsky et. al. 39 approximation 
are shown by blue dotted lines. 
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The coalescence radii p c were fitted for each com- 
posite particle in Ref. 75 to describe available data for 
the reaction Ne+U at 1.04 GeV/nucleon, but the fitted 
values turned out to be quite universal and were sub- 
sequently found to satisfactorily describe high-energy 
complex-particle production for a variety of reactions 
induced both by protons and nuclei at incident energies 
up to about 200 GeV/nucleon. These parameters used 
in the present version of CEM2k are: 

p c (d) = 90 MeV/c, 

Pc(t) = Pc( 3 He) = 108 MeV/c, (43) 
p c ( 4 He) = 115 MeV/c. 

If several cascade nucleons are chosen to coalesce into 
composite particles, they are removed from the sta- 
tus of nucleons and do not contribute further to such 
nucleon characteristics as spectra, multiplicities, etc. 

To describe fission-fragment production from heavy 
fissioning compound nuclei, we incorporated into the 
present version of CEM2k the thermodynamical fission 
model by Stepanov 77 with its own parameterizations 
for mass and charge widths, level-density parameters, 
fission barriers, etc. The already extended length of the 
current paper does not allow us to present here details 
of the thermodynamical fission model by Stepanov 77 
and sample results obtained with it. We plan to publish 
our fission results in a separate paper; for the moment, 
we only note that this model statistical by nature and 
coupled with CEM2k it provides similar results to those 
we got with the modified RAL model from GEM2, with 
some better agreement with the data for pre-actinides 
and some worse agreement for actinides. 

For brevity, we call the extended version of CEM2k 
described here as CEM2k2f (shorthand for "second 
fission model") and show below several results on 
complex-particle and fragment emission from several 
reactions, referring to them as from CEM2k2f. 

Fig. 19 shows examples of angle-integrated energy 
spectra of p, d, t, 3 He, and 4 He from 62 MeV pro- 
ton interactions with Al, 56 Fe, Y, and 120 Sn calculated 
with CEM2k2f (red histograms) compared with exper- 
imental data by Bertrand and Peelle 60 (symbols). For 
comparison, results obtained with previous versions of 
CEM are shown with black histograms, namely: cal- 
culations with CEM95 2 for 56 Fe, Y, and 120 Sn, and 
with CEM97 3 for Al. We still have some problems in a 
correct description of the high-energy tails of complex- 
particle spectra (that is probably related to direct pro- 
cesses of complex particle production like knock-out 
and pick-up which are not included in CEM2k2f), al- 
though the improved CEM2k2f provides a much bet- 
ter agreement with the data than the older versions 
CEM95 and CEM97. 

The new inverse cross sections and Coulomb barriers 
implemented here provide a correct and better descrip- 
tion of the evaporation parts of particle spectra, but 



require an adjustment of the "condensation" proba- 
bility 7j at the preequilibrium stage of reactions (see 
details in 1,68 ). CEM assumes 1 that during the pree- 
quilibrium part of a reaction pj excited particles (exci- 
tons) are able to condense with probability jj forming 
a complex particle which can be emitted during the 
preequilibrium state. In our first publication on CEM 1 , 
7j was roughly estimated as 

V-pjfa/A)*- 1 , (44) 

(A is the mass number of the pre-compound nucleus 
before emitting the particle j), and this estimation was 
used in all the succeeding versions of CEM. Generally, 
7j is a parameter of the preequilibrium model, and in 
the literature it is taken from fitting the theoretical 
preequilibrium spectra to the experimental ones, which 
gives rise to an additional, as compared to (44), de- 
pendence of the factor 7^ on pj, mass-number A, and 
excitation energy (see, e.g. Ref. 78 ). Our calculations 
show that one may obtain a reasonable description with 
CEM2k2f of all complex particles up to 4 He from many 
reactions if we use values for 7j - estimated using Eq. 
(44) divided by 3 for d and t, and by 4 for 3 He. The 
results shown in Fig. 19 are obtained with this estima- 
tion for 7j. It is clear that if we fit jj for each particle 
and each reaction as done in Ref. 78 , one may obtain a 
much better agreement with the data, but such exer- 
cises are outside the aim of our work. Our calculations 
of several reactions show that if we were to fit 7j - in 
CEM2k+GEM2, we could also obtain a much better 
description of complex particle spectra than the one 
shown in Fig. 14, but such a fit is as well beyond the 
aim of the present work. 

All spectra of particles from Al, 56 Fe, and 120 Sn 
shown on Fig. 19 are calculated without introducing 
a coalescence mechanism for complex-particle produc- 
tion. This is because for proton-induced reactions at 
these incident energies the multiplicity of fast nucle- 
ons emitted at the cascade stage of reactions is so 
small that the role of coalescence of complex parti- 
cles from such nucleons is negligible. As an example, 
in Fig. 19 we show only for Y (by blue histograms) 
calculated spectra that contain contributions from co- 
alescence production of complex particles. One can 
see that indeed the contribution from the coalescence 
mechanism is very small for this reaction and complex 
particle spectra calculated with and without this mech- 
anism almost coincide. 

Another example of the role of different mechanisms 
of particle production from these reactions is shown in 
Table 6, where we present the mean total multiplici- 
ties of complex particles from Y and 120 Sn calculated 
by CEM2k2f (last column) together with their contri- 
butions from preequilibrium, evaporation, and coales- 
cence reaction mechanisms (third to fifth columns). 
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Table 6. Role of different CEM2k2f reaction 
mechanisms in complex particle production 
(multiplicities) from p (62 MeV) + Y and 120 Sn 



Nucleus 


Ejectile 


Prec 


Evap 


Coales 


Total 




d 


2.96e-2 


1.41e-3 


1.83e-3 


3.28e-2 




t 


5.34e-3 


1.13e-4 


5.66e-6 


5.46e-3 




3 He 


1.92e-3 


1.13e-5 


2.83e-6 


1.94e-3 




4 He 


4.19e-2 


1.56e-2 




5.75e-2 


i2U Sn 


d 


3.17e-2 


6.51e-4 


1.75e-3 


3.41e-2 




t 


8.58e-3 


1.57e-4 




8.74e-3 




3 He 


7.59e-4 






7.59e-4 




4 He 


3.15e-2 


6.19e-3 




3.77e-2 



The percentage of the coalescence mechanism to the 
total complex particle yields is only 5.6% for d from Y 
and 5.1% for d from 120 Sn, while for t and 3 He from 
Y is only of the order of 0.1% and below the statistical 
errors of these reactions as simulated by CEM2k2f for 
4 He from Y and t, 3 He, and 4 He from 120 Sn. Nev- 
ertheless, our other studies show that at high incident 
energies 79 ' 80 and for heavy-ion induced reactions 75 , the 
coalescence mechanism of complex-particle production 
is important and often is the main mechanism of high- 
energy complex particle production; therefore we incor- 
porate this mode into CEM2k2f. 

It is informative to observe from Table 6 that due 
to high Coulomb barriers for complex particles from 
medium and heavy nuclei, the main contribution to 
production of such particles comes from preequilibrium 
emission, while the evaporation mode is less important 
for such reactions. 

Fig. 20 shows examples of 3 He and 4 He double- 
differential spectra from interaction of 300 and 480 
MeV protons with Ag calculated with CEM2k2f 
compared with experimental data by Green and 
Korteling 81 . One can see that CEM2k2f reproduces 
reasonably these data though it overestimates most of 
the spectra, especially the high-energy tails of them. 
Probably, by fitting jj for this reactions as discussed 
above one may obtain a much better agreement with 
the measurement, but we didn't do any additional fit- 
ting of 7j for these reactions. We also see some prob- 
lems with a correct description of angular distributions 
of 3 He and 4 He from these reactions and we plan to 
address this point in our further work. These results 
are preliminary and though we did not get a very good 
agreement with the data, they are encouraging to us 
as they help us to clarify ways of further improving 
CEM2k. 

From Fig. 20 we conclude that we probably get with 
CEM2k2f too much preequilibrium particle emission 
from these reactions. Such a conclusion may also be 
drawn from the results presented in Table 7, where we 
present a comparison of the total total yields of sev- 
eral particles and light fragments from 480 MeV p + 



Ag calculated with CEM2k2f with experimental data by 
Green, Korteling, and Jackson 82 . The results presented 
in Table 7 are preliminary, without any additional fit- 
ting of any parameters of CEM2k2f, nevertheless we 
see a reasonable agreement with the measurement not 
only for production of 3 He and 4 He, but also for heavier 
fragments. We calculated this reaction with CEM2k2f 
both including preequilibrium-particle emission (third 
column in Table 7) and ignoring it (last column in Ta- 
ble 7). We can see that for many of the fragment- 
production cross sections ( 8 Li, 9 Li, 7 Be, 9 Be, and n Be 
shown in the table) the experimental yields are higher 
than the ones calculated without preequilibrium emis- 
sion, and lower than those with preequilibrium emis- 
sion. Though these results are preliminary and will 
change when we develop further our model, they may 
serve as a rough indication that in the present version 
of CEM2k2f we have too much preequilibrium parti- 
cle emission, just as we got with CEM2k+GEM2 dis- 
cussed in Sec. III. We recall again that in 56 we address 
this point both in our CEM2k 10 and LAQGSM 79 codes 
merged with GEM2 37 ' 38 using an approach based on 
Ref. 67 , and the results obtained there are promising. 

Table 7. Comparison of the measured 82 total pro- 
duction cross sections (mb) for several ejectiles from 
480 MeV p + Ag with our preliminary CEM2k2f re- 
sults with and without preequilibrium emission 



Ejectile 


Exp. data 


CEM2k2f 


CEM2k2f 






with Prec 


no Prec 


n 




3.97e+3 


4.50e+3 


P 




9.30e+2 


9.94e+2 


d 




2.85e+2 


3.51e+2 


t 




1.27e+2 


1.35e+2 


3 He 


3.42e+l 


3.98e+l 


3.02e+l 


4 He 


4.44e+2 


7.68e+2 


l.lle+3 


6 He 




6.60e+0 


3.61e+0 


8 He 




1.31e-l 


1.61e-2 


6 Li 


3.67e+l 


9.91e-2 


9.28e+0 


7 Li 


4.20e+l 


3.59e-2 


5.27e+0 


8 Li 


4.48e-l 


1.24e-3 


5.31e-l 


9 Li 


6.98e-2 


1.08e-4 


8.37e-2 


7 Be 


8.36e-l 


5.85e-3 


1.20e+0 


9 Be 


9.22e-l 


1.06e-2 


1.73e+0 


10 Be 


4.16e-l 


6.23e-4 


1.52e-l 


n Be 


1.29e-2 


2.24e-5 


1.66e-2 


12 Be 




7.65e-7 


1.02e-3 


8 B 




5.98e-5 


3.46e-2 


10 B 


4.85e-l 


1.35e-3 


4.22e-l 
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Figure 19. Angle-integrated energy spectra of p, d, t, 3 He, and 4 He from 62 MeV protons on Al, 56 Fe, Y, and 
120 Sn. Calculations from the CEM2k2f code are shown by red histograms and from CEM95 2 (for 56 Fe, Y, and 
120 Sn) and CEM97 3 (for Al) are shown by black histograms. For Y, the blue histograms show CEM2k2f results 
taking into account contributions from coalescence of complex particles from fast cascade nucleons, as described 
in the text. Experimental data (symbols) are by Bertrand and Peelle 60 . 
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Figure 20. Double-differential spectra of 4 He and 3 He from 480 and 300 MeV protons on Ag. Calculations from 
the CEM2k2f code are shown by solid histograms and contribution from only the evaporation mode, by dashed 
histograms. Experimental data (symbols) are by Green and Korteling 81 . 
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Finally, we have investigated the well known code 
GEMINI by Charity 83 as an alternative way to describe 
production of various fragments by merging GEMINI 
with both our CEM2k and LAQGSM. The preliminary 
results we find for spallation, fission, and fragmentation 
products from several reactions are promising and we 
will present our results from this study in a separate 
paper. 

V. FUTURE WORK 

The two modifications of the CEM2k model presented 
here allow us to describe satisfactorily many fission 
and fragmentation reactions in addition to the spal- 
lation reactions which are already quite well described 
by CEM2k. 

Our results presented in Sec. Ill show that the 
Generalized Evaporation Model code GEM2 by Furi- 
hata is a useful tool and when merged with CEM2k 
(or LAQGSM) it allows us to reproduce a large vari- 
ety of fission and fragmentation measurements, though 
some tuning of several GEM2 parameters will be re- 
quired. We may choose a model similar to GEM2 in 
our codes to describe fission, fragmentation, and evap- 
oration yields, but it must be significantly extended and 
improved in order to properly describe complex parti- 
cles, light and heavy fragments not produced in fission, 
and more importantly, to become self-consistent and 
better grounded from a physical point of view. 

The latest version of CEM2k described in Sec. IV, 
CEM2k2f, contains new inverse cross sections and 
Coulomb barriers, new routines to calculate widths of 
emitted particles and to simulate their kinetic energies 
using arbitrary approximations for the inverse cross sec- 
tions and Coulomb barriers, allows evaporation of up 
to 66 different particles both in the ground and excited 
states and emission of up to 29 different preequilib- 
rium particles, takes into account coalescence of com- 
plex particles from fast nucleons emitted at the cas- 
cade stage, and describes much better complex-particle 
spectra and yields than its precursors. 

Our work is not finished yet. For example, at the 
evaporation and preequilibrium stages we still need to 
make an "ad hoc" decision about how many different 
particles should be included into consideration so that 
we have production of light fragments heavier than 4 He 
but maintain the running time of the code at a reason- 
able level. Our new algorithms and routines do not 
impose any limitations to the energy dependence of 
the particle emission probabilities, and we plan to try 
to improve the present description of the level densities 
and particle emission widths with a hope to describe 
better the tails of complex-particle spectra. 

We consider also implementation into CEM2k a 
model of fragmentation, to be able to describe emission 
of fragments like Na, Mg, and heavier, and possibly, 
incorporation of the Fermi break-up model to replace 



the preequilibrium and evaporation models when the 
residual excited nucleus after the cascade stage of a 
reaction is very light. 

Our work on fission in CEM2k is also not finished. 
For instance, we are not satisfied with the situation 
that in both the improved versions of the CEM2k dis- 
cussed here we still have an additional input parameter 
to describe fission cross sections: either C(Z), in the 
GEM2/RAL approach for nuclei with Z > 88, and/or 
af/a n , in the same model for lighter fissioning nuclei 
with Z < 88 and when using the model by Stepanov. 
In addition, to get with CEM2k a better description of 
fission-fragment distributions, the approximations used 
in both GEM/RAL and Stepanov's fission models to 
simulate A, Z, and kinetic-energy distributions for fis- 
sion fragments should be further adjusted in our code. 

We emphasize that it is not sufficient to analyze 
only A and Z distributions of the product yields when 
evaluating the type of model discussed here, as is often 
done in the literature, but that it is essential to study 
all the separate isotope yields as well as the spectra of 
light particles and fragments. This should be done for 
different target nuclei, different incident particles, and 
at different energies. 
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